### Set working directory
setwd("P:/Fabian/PAPER_Popgen_SouthAfrica/Submission")

### Load libraries
library(ggplot2)
library(MASS)
library(rmarkdown)
library(knitr)
library(lme4)
library(MuMIn)
library(ggthemes) # Load
library(datasets)
library(plyr)
library(dplyr)
library(tidyverse)
library(coin)

### Analysis of all loci on all samples
### Linkage disequilibrium: BY correction for all loci 
p.adjust(c(0.319,	0.5081,	1,	0.0007,	0.9993,	0.6246,	0.8635,	0.7747,	0.8122,	0.873, 0.2959,	0.6512,	0.1129,	0.2435,	0.009,	0.0697,	0.0105,	
           0.0821,	0.9657,	0.2792,	0.202,	0.0223,	0.5546,	0.5939,	0.0197,	0.3961,	0.2169,	1,	0.3314,	0.9622,	0.383,	0.1596,	0.4174,	0.2695,
           0.6003,	1,	0.4111,	0.595,	0.2757,	0.3482,	0.2584,	0.5225,	0.4197,	0.8574,	0.9827,	0.065,	0.4618,	0.5655,	0.8585,	0.3999,	0.3344,	0.857,
           0.4722,	0.8832,	0.6736), method="BY")


### Correlation between FIS and FST for all loci on all samples
FISFST <- read.csv("FISFST_allloci_allsamples.csv", sep=",", row.names=NULL)
FISFST
with(FISFST, cor.test(FST, FIS, alternative="greater", method="spearman"))


### Correlation between FIS and number of blanks for all loci on all samples
FISblanksall <- read.csv("FISblanks_allloci_allsamples.csv", sep=",", row.names=NULL)
FISblanksall
with(FISblanksall, cor.test(Nb, FIS, alternative="greater", method="spearman")) 


### Correlation between FIS and number of blanks for 7 loci (Gb5, 28, 35, 66, 73, 92, 165) on all samples
FISblanks7lociall <- read.csv("FISblanks_7loci_allsamples.csv", sep=",", row.names=NULL)
FISblanks7lociall
with(FISblanks7lociall, cor.test(Nb, FIS, alternative="greater", method="spearman")) 

### Correlation between FIS and number of blanks for 4 loci (Gb 5, 35, 66, 73) on all samples
FISblanks4lociall <- read.csv("FISblanks_4loci_allsamples.csv", sep=",", row.names=NULL)
FISblanks4lociall
with(FISblanks4lociall, cor.test(Nb, FIS, alternative="greater", method="spearman")) 


### Short allele dominance 
SAD <- read.csv("SAD.csv",sep=",", row.names=NULL)
SAD

SAD_Gb28 <- subset(SAD,loci=="Gb28")
SAD_Gb28
with(SAD_Gb28, cor.test(Allele, Capf, alternative="less", method="spearman"))


SAD_Gb35 <- subset(SAD,loci=="Gb35")
SAD_Gb35
with(SAD_Gb35, cor.test(Allele, Capf, alternative="less", method="spearman"))

SAD_Gb48 <- subset(SAD,loci=="Gb48")
SAD_Gb48
with(SAD_Gb48, cor.test(Allele, Capf, alternative="less", method="spearman"))
with(SAD_Gb48, cor.test(Smallf, Weights, alternative="less", method="spearman"))

lm_SAD_Gb48 <- lm(Smallf ~ Allele, data=SAD_Gb48, Weights=weight)
summary(lm_SAD_Gb48) 

SAD_Gb5 <- subset(SAD,loci=="Gb5")
SAD_Gb5
with(SAD_Gb5, cor.test(Allele, Capf, alternative="less", method="spearman"))

SAD_Gb66 <- subset(SAD,loci=="Gb66")
SAD_Gb66
with(SAD_Gb66, cor.test(Allele, Capf, alternative="less", method="spearman"))

SAD_Gb70 <- subset(SAD,loci=="Gb70")
SAD_Gb70
with(SAD_Gb66, cor.test(Allele, Capf, alternative="less", method="spearman"))

SAD_Gb72 <- subset(SAD,loci=="Gb72")
SAD_Gb72
with(SAD_Gb72, cor.test(Allele, Capf, alternative="less", method="spearman"))

SAD_Gb158 <- subset(SAD,loci=="Gb158")
SAD_Gb158
with(SAD_Gb158, cor.test(Allele, Capf, alternative="less", method="spearman"))

SAD_Gb73 <- subset(SAD,loci=="Gb73")
SAD_Gb73
with(SAD_Gb73, cor.test(Allele, Capf, alternative="less", method="spearman"))

SAD_Gb92 <- subset(SAD,loci=="Gb92")
SAD_Gb92
with(SAD_Gb92, cor.test(Allele, Capf, alternative="less", method="spearman"))

SAD_Gb165 <- subset(SAD,loci=="Gb165")
SAD_Gb165
with(SAD_Gb165, cor.test(Allele, Capf, alternative="less", method="spearman"))

### Linkage disequilibrium: BY correction for 9 selected loci
p.adjust(c(0.3265,	0.5193,	0.0006,	0.6329,	0.864,	0.7794,	0.821,	0.8735,	0.2963,	0.1097,	0.0085,	0.0715,	0.0115,
           0.081,	0.9701,	0.2031,	0.5426,	0.5946,	0.0201,	0.4044,	0.2177,	1,	0.3921,	0.5911,	0.269,	0.3584,	0.0654,
           0.4602,	0.5638,	0.8668,	0.4043,	0.3298,	0.8541,	0.4636,	0.8868,	0.6678), method="BY")

### Correlation between FIS and FST for 9 selected loci on all samples
FISFST9loci <- read.csv("FISFST_9loci_allsamples.csv", sep=",", row.names=NULL)
FISFST9loci
with(FISFST9loci, cor.test(FST, FIS, alternative="greater", method="spearman"))

### Correlation between FIS and FST for 9 selected loci on field samples only
FISFST9locifieldsamples <- read.csv("FISFST_9loci_fieldsamples.csv", sep=",", row.names=NULL)
FISFST9locifieldsamples
with(FISFST9locifieldsamples, cor.test(FST, FIS, alternative="greater", method="spearman"))

### Correlation between FIS and number of blanks for selected 9 loci on field samples
FISblanks9loci_fieldsamples <- read.csv("FISblanks_9loci_fieldsamples.csv", sep=",", row.names=NULL)
FISblanks9loci_fieldsamples
with(FISblanks9loci_fieldsamples, cor.test(Nb, FIS, alternative="greater", method="spearman")) 

### Correlation between FIS and number of blanks for loci with null alleles (Gb5, 35, 66, 73, 92)
FISblankslocinonulls_fieldsamples <- read.csv("FISblanks_locinonulls_fieldsamples.csv", sep=",", row.names=NULL)
FISblankslocinonulls_fieldsamples
with(FISblankslocinonulls_fieldsamples, cor.test(Nb, FIS_nulls, alternative="greater", method="spearman")) 


